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As it is known from visible light experiments, silicon under femtosecond pulse irradiation can 
undergo the so-called ’nonthermal melting’ if the density of electrons excited from the valence to 
the conduction band overcomes a certain critical value. Such ultrafast transition is induced by 
strong changes in the atomic potential energy surface, which trigger atomic relocation. However, 
heating of a material due to the electron-phonon coupling can also lead to a phase transition, 
called ’thermal melting’. This thermal melting can occur even if the excited-electron density is 
much too low to induce non-thermal effects. To study phase transitions, and in particular, the 
interplay of the thermal and nonthermal effects in silicon under a femtosecond x-ray irradiation, we 
propose their unified treatment by going beyond the Born-Oppenheimer approximation within our 
hybrid model based on tight binding molecular dynamics. With our extended model we identify 
damage thresholds for various phase transitions in irradiated silicon. We show that electron-phonon 
coupling triggers the phase transition of solid silicon into a low-density liquid phase if the energy 
deposited into the sample is above ~ 0.65 eV per atom. For the deposited doses of over ~ 0.9 eV 
per atom, solid silicon undergoes a phase transition into high-density liquid phase triggered by an 
interplay between electron-phonon heating and nonthermal effects. These thresholds are much lower 
than those predicted with the Born-Oppenheimer approximation (~ 2.1 eV/atom), and indicate a 
significant contribution of electron-phonon coupling to the relaxation of the laser-excited silicon. We 
expect that these results will stimulate dedicated experimental studies, unveiling in detail various 
paths of structural relaxation within laser-irradiated silicon. 

PACS numbers: 41.60.0, 64.70.K-, 42.65.Re, 61.80.Ba 


I. INTRODUCTION 

Nonthermal melting is a well-established concept, known for over two decades both theoretically M and ex¬ 
perimentally Thus, it can be surprising how well ’thermal’ models could sometimes reproduce experimental 

observations of phase transitions triggered by laser pulse irradiation of solids fl0Ml~3l ] . This controversy stimulated 
intense discussions (Til, [71 . It originates from the fact that the two approaches: thermal (relying on electron-phonon 
heating of the atomic system on an unchanged potential energy surface) and nonthermal (describing changes of in¬ 
teratomic potential surface via electron excitation while excluding electron-phonon coupling), are based on different 
assumptions and approximations, which are rarely studied together due to prohibitive computational complexity. 

The nonthermal effects in solids are typically studied within the Born-Oppenheimer approximation [Ij-Q. The 
models treating electron-phonon coupling, such as atomistic-continuum models, TTM-MD E2El,E!Ezl include the 
non-adiabatic effects by adding an empirical electron-phonon coupling parameter. To our knowledge, there were only 
a few attempts to incorporate both effects for solids. For example, such attempts were made in a phenomenological 
manner in Refs. El El- 

It is generally believed that for a low deposited dose, thermal melting of a semiconductor or an insulator can be 
induced, while for a higher dose, when typically ~10% of the valence-band electrons are excited to the conduction 
band 0E1, a nonthermal melting occurs. This is, however, not always the case: diamond is a counterexam ple. For 
diamond, at a lower deposited dose, the nonthermal graphitization occurs and not the thermal amorphization |19j, [20j. 
Thus, one cannot say a priori which mechanism dominates for a specific material, and a dedicated analysis is required 
for each case. 

Apart from the conventional visible-li ght lasers, the 4 th generation light sources, the free electron lasers (FEL, 
such as FLASH [U, LCLS [U, SACLA [23|, FERMI [H]), emitting intense femtosecond x-ray pulses can shed new 
light on the problem. Since almost a decade, FELs have stimulated rapid advances in many scientific fields. Damage 
of semiconductors under a femtosecond x-ray irradiation starts with photoabsorption, which excites electrons from 
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the valence band or, in contrast to a visible-light irradiation, from deep atomic shells, to high-energy states of the 
conduction band 1251 1271 . The deep-shell holes (K- and L-shell holes for silicon) then decay via Auger processes. This 
is the dominant relaxation channel for light elements (28]. Auger decay of a hole leads to an excitation of an electron 
from a higher shell to the conduction band, with the transfer of the excess energy from a relaxing deep-shell hole to 
the excited electron. The ejected photo- and Auger-electrons scatter further via inelastic channels (impact ionization 
of valence band or deep-shell electrons), or elastic channels (scattering on atoms or phonons). The impact-ionization 
cascading is the most important relaxation process for high-energy electrons. It typically occurs on a femtosecond 
timescale, and finishes when the electron loses its energy below the impact ionization threshold [2?j [30]. In contrast, 
the elastic electron-phonon scattering dominates for low-energy electrons, leading to significant electron energy losses 
only at longer (typically picosecond) timescales [13, UK, 32:]. Apart from that, the transiently excited state of the 
electronic subsystem induces a change of the atomic potential energy surface. This can lead to nonthermal melting 
described above, i.e., in covalently bonded semiconductors, the enforced population of antibonding states within the 
conduction band can trigger an ultrafast rearrangement of atoms which then attempt to minimize the potential energy. 

To study phase transitions in silicon triggered by a femtosecond FEL pulse, we use our recently developed hybrid 
model [19L l20j ] which traces nonequilibrium kinetics of electrons under ultrashort laser irradiation and the following 
rearrangement of atoms. The model relies on: (i) a Monte Carlo description of nonequilibrium high-energy elec¬ 
trons, (ii) Boltzmann kinetic approach for low-energy electrons, and (iii) tight-binding molecular dynamics (TBMD) 
method tracing atomic trajectories, transient electronic band structure, and evolution of the interatomic potential 
energy surface. The original TBMD method proposed by Jeschke et al. [33, HU was based on the Born-Oppenheimer 
approximation. However, to address thermal melting occurring via electron-lattice coupling, it is necessary to extend 
the model beyond the Born-Oppenheimer approximation j3f|. We accordingly modify our hybrid model [l9l |20| by 
including the non-adiabatic coupling between electrons and the lattice, which is usually known as ’electron-phonon 
coupling’. The proposed model calculates at each time step the respective transition rate from the matrix element for 
electron-atom (ion) coupling known in ab-initio femto-chemistry |36l l37|. The rate is then included as a Boltzmann 
collision integral |27i l32j in the evolution equation for electron distribution. The extended model thus accounts for 
both thermal and nonthermal effects, allowing us to study their occurrence within one consistent theoretical frame¬ 
work. Note also that the same approach can be used in any ab-initio model, such as, e.g., density-functional-theory 
molecular dynamics. 


II. MODEL 
A. Hybrid model 

The hybrid model addressing processes occurring in a semiconductor during its irradiation with VUV-rays or x-rays 
combines various simulation methods. It was developed in [jj| and described in detail in Ref. [38]. The basic ideas 
are as follows. The atom dynamics is traced with the classical molecular dynamics simulation method (MD) with 
periodic boundary conditions [l9l. [ikil. [39|. This method requires a knowledge of the potential energy surface, which 
determines the forces acting on the atoms. 

The potential energy surface together with the transient electronic band structure are calculated by a direct diago- 
nalization of a tight binding (TB) Hamiltonian. The Hamiltonian is following the evolution of the atomic configuration 
within a simulation box, thus, changing in time. The forces acting on atoms and the electron-lattice (electron-phonon) 
energy exchange depend additionally on th e sp ecific state of the electronic subsystem. The corresponding potential 
energy surface is calculated as in Refs. |33ll39|: 

= fe(Ej, t)Ej + -Erepdrjj}) , (1) 

i 

where the repulsive part is describing the effective repulsion of atomic cores, E rep ({r,;j}). This potential energy is 
used within the Parrinello-Rahman Lagrangian under the constant-pressure condition. The corresponding equations 
of motion are presented in [u| [H, j4lj • 

The transient electron distribution function f e (Ei,t) enters Eq. Jl]). Therefore, one has to trace simultaneously 
the evolution of the state of electronic ensemble. Additionally, the transient energy levels Ej are obtained by a 
diagonalization of the tight binding Hamiltonian (for more details see 0 )- 

After a VUV- or x-ray irradiation, the transient electron distribution function has a shape of the so-called ’bump- 
on-hot-tail’-distribution [25l. 12(1 Pill Phil . This typical shape consists of nonthermalized high-energy electrons, and of 
(nearly-) thermalized low-energy electrons within the valence and the bottom of the conduction band. Utilizing this 
fact, we split the electron ensemble in two parts, treating each of them with a dedicated (computationally efficient) 


3 


method. The Monte Carlo (MC) method is used to describe the transient nonequilibrium kinetics of high-energy 
electrons and their secondary cascading as well as the photoabsorption and Auger decays of atomic deep-shell holes 
[ISilSIilil. More details on the Monte Carlo model and the cross sections used can be found in [HI iH| . The 
cross sections used for electron scattering in silicon can be found in Ref. [44j. 

A simplified Boltzmann equation is applied to describe low-energy electrons. The high-energy-electron and the 
low-energy-electron domains are interconnected, as electrons can gain or lose energy and go from one domain to 
another. This forms the source/sink terms for the low-energy part [45l . 146} , as the changing number and energy of 
low-energy electrons affect directly their distribution. Additionally, atomic motion and the evolution of the electronic 
band structure also influence the temperature of electrons and their chemical potential fl9t |. 

We developed a dedicated technique to treat electron-atom energy exchange via non-adiabatic channel (electron- 
phonon coupling), as it will be explained in detail in section III Bl Such combined treatment enables us to trace 
modification of the atomic potential caused by the excitations of electrons and the electron-atom energy exchange, 
addressing possible thermal and nonthermal phase transitions simultaneously. 


B. Boltzmann equation for low-energy electrons 


Evolution of the electron distribution function on the energy levels obtained from the diagonalization of the TB- 
Hamiltonian can be traced by means of Boltzmann collision integrals [47} : 


dh 

dt 




+ s, 


3 3 


( 2 ) 


where 7?“ e is an electron-electron collision integral, is an electron-atom collision integral, and S' is a source term, 

describing the electrons arrivingand leaving the low-energy domain, in accordance with the introduced separation of 
low- and high-energy domains [l9j, [45], [48| . 

As the low-energy electrons are in a nearly thermalized state already after a few femtoseconds since the beginning of 
the laser pulse, the electron-electron collision integral turns to zero. In our model, we ensure the electron distribution 
function to be a Fermi-function due to the assumed instant thermalization of electrons at each time-step, similarly to 
Ref. jl9} . Note that any possible slight deviation of the exact distribution function from the equilibrium Fermi-shape 
affects only negligibly the atomic motion and the phase transition [49}. 

For the Fermi distributions one can define the corresponding chemical potential and electronic temperature 0. 
This temperature is generally different from the atomic temperature, thus, electron-atom scattering integral governs 
the energy flow, Q, between the two systems: 




—at 
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(3) 


The electron-atom scattering integral depends on the transient distributions of electrons and of atoms, and on the 
matrix element describing their interaction. Within the Born-Oppenheimer approximation, an atomic motion does 
not trigger any electron transition between the energy levels, as electrons are assumed to adjuste instantly to a new 
configuration. Thus, the approximation is inherently incapable of reproducing electron-atom energy exchange, and 
the collision integral I^~ at = 0. In order to trace electron-atom coupling, non-adiabatic effects (electron-phonon 
coupling) must be explicitly included. Generally, electron-atom collision integral can be written in the form (similar 

to [53): 


N 9-tt N 

= TE \ M e-at( E i, E j)\ 2 X (4) 

3 = 1 3 =1 

7e(£i)(2 - fe(Ej)) - /e(F,)(2 - f e ( E i))gat( E i - E j) , 
for i > j, 

' f e (Ei )(2 - f e (E j ))g at (E j - Ei) - / e (F,)(2 - / e (F,)) , 

„ for i < j, 

where g at is the integral of the atomic distribution function, and M e _ at (E i: Ej) is the matrix element for electron-atom 
(ion) scattering. We calculate the matrix element with the method used in non-adiabatic molecular dynamics applied 
in femtochemistry [§6[ [37]: 

Af e -at ( E i i E j) = 
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FIG. 1: Electron-lattice thermalization: comparison of the non-adiabatic scheme with the Born-Oppenheimer approximation. 
Initial system conditions are: the electron temperature T e = 10000 K, the atomic temperature T a t = 300 K. A super-cell of a 
constant volume with 216 atoms is used. 


\ «*(* - St )\j(t)) - 0'(* ~ St)\i{t))) (Ej - Ei ) , (5) 

where M e - at (Ei , Ej ) is the matrix element for electron transition between the levels Ei and Ej induced by the atomic 
motion, taken as a mean value on the current and previous time-steps: Ei = ( Ei(t) + Ei(t — St ))/2 ; | i(t)) is the 
electron wave-function at the time-instance t obtained as an eigenfunction of the TB-Hamiltonian. Note that the 
electron wave-functions from two sequential time-steps (t — St and t) are entering Eq. (0 [35j. The derivation of the 
collision integral Eq.fl]), the matrix element, Eq.flSJ, and the numerical details of their calculations are presented in 
Appendix. 

Once we calculate the average heat flow between the electrons and ions, Eq.®, we use the velocity scaling for atoms 
to accommodate the excess energy transferred from (or to) electrons on the current timestep [HU • The corresponding 
change of the electron energy is introduced to the electron distribution function. 

A test-case is presented in Fig[T| It shows the relaxation of the electron-atom system at an initial nonequilibrium 
between hot electrons and room-temperature atoms. Electron and atom temperatures are T e = 10000 K and T a = 300 
K, correspondingly. The number of atoms within a unit cell is 216. It is sufficiently large so that there is no 
artificial influence of the number of atoms on the thermalization timescale (detailed convergence study is presented 
in Appendix). The Born-Oppenheimer approximation does not allow for any energy exchange between the electrons 
and lattice, while the non-adiabatic scheme yields reasonable timescales for electron-lattice thermalization, similar to 
the empirical estimates Ennm. 


III. RESULTS 

A. Interplay of thermal and nonthermal effects 

With the extended model we have simulated evolution of laser-irradiated silicon at various radiation doses absorbed 
per atom. 

FigH] shows snapshots of the atomic positions of silicon after an FEL pulse of 10 fs duration, Hoj =1 keV, and the 
absorbed dose of 0.7 eV/atom. At such deposited doses, silicon reaches only the low-density liquid phase (LDL) 
[HHj, characterized by an electronic phase transition into a semi-metallic state with closed band gap (below). Atomic 
temperature then exceeds silicon melting temperature of 1687 K (FigJS]). The local order in the atomic structure is 
preserved. Silicon remains in the LDL state during the whole simulation time, i. e., up to 50 ps (not shown). 

The snapshots of atomic positions within silicon irradiated with an FEL pulse of 10 fs duration, huj= 1 keV, and the 
absorbed dose of 0.9 eV/atom are shown in Figj3] After absorbing this dose, silicon reaches the high-density liquid 
(HDL) phase with amorphization f53| . This is a result of the interplay between thermal heating and nonthermal 
changes in the interatomic potential. They trigger atomic relocations on an short time-scales of 300 — 500 fs. These 
time-scales match very well the experimental observations of nonthermal melting sum. Similarly to observed 
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FIG. 2: Transition to LDL phase: snapshots of atomic positions in silicon irradiated with 10 fs laser pulse of hw=l keV photon 
energy at the absorbed dose of 0.7 eV/atom: (a) t = 0 fs, (b) t = 300 fs, (c) t = 0.5 ps, and (d) t = 1 ps. X, Y, and Z axes 
are shown (left-bottom of each panel). 



FIG. 3: Transition to HDL phase: snapshots of atomic positions in silicon irradiated with 10 fs laser pulse of ftw =1 keV photon 
energy at the absorbed dose of 0.9 eV/atom: (a) t = 0 fs, (b) t = 300 fs, (c) t = 0.5 ps, and (d) t = 1 ps. X, Y, and Z axes 
are shown (left-bottom of each panel). 


in [H, the HDL phase is reached after an intermediate LDL phase. 

Figure @] shows how the volume of the Parrincllo-Rahman super-cell changes after the irradiation with different 
fluences. As the number of atoms and ions within the super-cell is conserved, the decrease or increase of the super¬ 
cell volume corresponds to the increase or decrease of the atomic density respectively. For the doses of 0.5 and 0.7 
eV per atom which are below the nonthermal damage threshold, one can see an increase of the super-cell volume, 
corresponding to the decrease of the atomic density only during the electron cascading time (~first hundreds of fs). 
The dose of 0.7 eV/atom corresponds to the phase transition to the LDL phase, as indicated by the band gap collapse 
in FiglU The doses above the nonthermal melting threshold, 0.9 and 1 eV/atom in FigQ] induce a transient increase 
of the volume with its shrinkage later. This reflects the intermediate transition of silicon to the LDL phase followed 
by the transition to HDL phase. 

Fig.[5]presents the band gap of silicon after different energy depositions. Band gap is defined as the energy difference 
between the closest eigenstates above and below the Fermi energy. Band gap width shrinks to nearly zero already at 


6 



FIG. 4: Volume of the super-cell of silicon (216 atoms) irradiated with 10 fs laser pulse of hu )=1 keV photon energy at various 
absorbed doses. 



Time (fs) 


FIG. 5: Band gap of silicon irradiated with 10 fs laser pulse of Tiw= 1 keV photon energy plotted as a function of time at various 
absorbed doses. 


~ 0.7 eV/atom, showing that silicon is in a semi-metallic state (LDL). This is in agreement with our results indicating 
that the phase transition into LDL phase occurs above the threshold of ~ 0.65 eV/atom (Fig. [2]). 

Fig. [6] shows that deep-shell holes decay quickly via Auger-decay, as discussed above, and their energy is brought 
back to the electronic system on a sub-100 fs scale. Electron cascading finishes within ~ 250 — 300 fs (also seen in 
Fig.©, and the energy is transferred to the low-energy electrons. The figure [7] confirms that the energy redistribution 
between different subsystems of the irradiated material occurs on femtosecond scale, starting with the excitation of 
electrons and holes. The total energy is conserved, as we did not consider any energy transport from the system, 
assuming periodic boundaries. In contrast to the diamond-to-graphite phase transition reported in HHH , the atomic 
potential energy does not exhibit a rapid jump at the beginning of phase transition. Instead, the system is relaxing on 
picosecond timescales. For higher absorbed doses (0.9 and 1 eV/atom in figure [3), silicon turns into the high-density 
liquid phase, which is also reflected by the potential energy curve: it is slightly raising (inset in the Fig. Q, while 
atomic temperature is decreasing (see below in Fig. ©. The timescales of these changes match the timescales of the 
super-cell volume contraction shown in Fig. [4] 

Atomic and electronic temperatures for the same absorbed doses are shown in Fig. [5] The electron temperature 
is increasing during the laser pulse and during the relaxation of high-energy electrons and deep-shell holes (first 200 
fs). Later, while electrons transfer their energy to the lattice, the electron temperature is decreasing. This takes a 
few picoseconds. After that time, electrons are equilibrated with the atoms. The electron temperature oscillations 
(especially pronounced for low dose irradiation) are caused by the super-cell volume oscillations, hence, they can be 




















Time (fs) 

FIG. 6: Densities of high-energy electrons (top panel) and holes in L\ and 1/2,3-shell holes (bottom-panel) as a function of time. 
They are expressed as a percentage of the initial valence-electron density. 


considered as a simulation artifact. Such oscillations are reflected either in electron temperature (for fixed energy), 
or in the energy (for fixed temperature); we chose the first scheme in the simulation. Atomic temperature oscillations 
reflect physical process: the exchange between kinetic and potential energies of atoms. At the highest fluence, 
atomic temperature increases rapidly already within ~ 300 fs. This reflects a strong interplay between thermal and 
nonthermal effects within the system. 

The number of low-energy conduction-band electrons is shown in Fig. [9] For the absorbed doses considered above 
it never reaches the critical value of 9% which purely nonthermal models predict [I, 11 • However, the number of 
electrons is sufficiently high to trigger a phase transition, which is then due to the thermal heating of the system with 
nonthermally-weakened interatomic bonds. To mention, the peak conduction-band electron densities are close to the 
ones estimated experimentally Q. 

Note that the density of the excited electrons reached after FEL irradiation at the pulse fluences considered here, 
which provide absorbed dose on the level of ~ 1 eV/atom, is of the order of a few percent of the solid density 
(10 21 — 10 22 cm -3 )^ This is a high density for an electron plasma. Its (partial) thermalization is then known to be 
very rapid I25l 12(1 I41 h 43| . confirming the assumptions made above. On the other hand, being only a few percent of 
the solid density, it ensures that the applicability condition of the tight binding scheme (low excitation regime) is not 
violated. 


B. Damage threshold for silicon as a function of photon energy 

Following the procedure from Ref. [20], we estimate the damage threshold of silicon for different photon energies. 
We checked (not shown) that the damage threshold in terms of deposited energy per atom is almost independent of the 
incoming photon energy. Thus, by converting the dose per atom into the units of incoming fluence with one-photon 
absorption cross section from Refs. [55l457j . we obtain the damage threshold predictions shown in Fig. [10] They can 
be directly verified experimentally. No effects of thermal diffusion and particle diffusion were taken into account for 
the predictions from Fig. [TO] They can play a role for the case of small skin-depth (the minimum around 20 eV on 
the picture). Due to these effects, and possible re-solidification governed by the energy flows out of the laser spot, at 
such photon energies our calculations might underestimate the experimentally measured damage thresholds. 


C. Purely nonthermal melting of silicon 

Finally, for comparison we show the predictions for silicon melting obtained after excluding the non-adiabatic effects 
(electron-phonon coupling). Within the adiabatic (Born-Oppenheimer) scheme the calculated damage threshold for 
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FIG. 7: Energy redistribution in silicon irradiated with 10 fs FEL pulse of fku= 1 keV photon energy at the absorbed doses of 
0.5, 0.7, 0.9 and 1 eV/atom. The black solid line is the potential energy of atoms; red dashed line is the total energy of atoms 
(potential and kinetic); blue dotted line is the total energy of atoms and electrons (energy of the system excluding deep-shell 
holes); and the green dash-dotted line is the total energy of the system. The inset is showing longer timescales behavior of the 
potential energy of atoms for the case of 0.9 eV/atom absorbed dose. 


nonthermal melting of silicon appears to be at the absorbed dose of 2.1 eV/atom. This corresponds to 9% of electrons 
excited from the valence to the antibonding states of the conduction band. This threshold value of the electron density 
is in an excellent agreement with earlier works , based on the adiabatic scheme. 

The atomic snapshots shown in Fig. [Til demonstrate the final state to be an amorphous high-density liquid [58| . 
The transition to that state proceeds on sub-picosecond timescales at which thermal effects - if included - would play 
a role. The intermediate LDL phase can be identified by the collapsed band gap and increased volume of the modeled 
super-cell (decreased density) prior to the final contraction to the HDL phase. 

As discussed above, the model of nonthermal melting of silicon based on the Born-Oppenheimer approximation does 
not allow for electron-phonon coupling, i.e., the excited electrons do not exchange energy with atoms. The damage 
then occurs as a purely nonthermal effect of weakening the interatomic bonds due to electron excitation. This results 
in much higher damage threshold than that obtained with non-adiabatic approach. These finding also indicates that 
thermal models such as EMMS can be accurately applied for silicon irradiated with laser pulses of low fluences, 
when the nonthermal effects do not play a significant role. 


IV. CONCLUSION 

In the present work, we studied phase transitions in silicon under a femtosecond irradiation with x-ray laser. In 
order to account both for thermally and nonthermally triggered transitions, we have extended our recently developed 
hybrid model fl9l . [2l| by including non-adiabatic electron-phonon coupling. In this way heating of a material due to 
the electron-phonon coupling can also be treated. The developed scheme is general and can be used in any ab-initio 
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FIG. 8: Electronic and atomic temperatures in silicon irradiated with 10 fs FEL pulse of fku=l keV photon energy at the 
absorbed doses of 0.5, 0.7, 0.9 and 1 eV/atom. Silicon melting temperature is 1687 K. 



FIG. 9: Density of low-energy conduction-band electrons in silicon after an FEL pulse of 10 fs duration, hw=l keV photon 
energy at the absorbed doses of 0.5, 0.7, 0.9 and 1 eV/atom. 


molecular dynamics model. 

We demonstrated that for silicon under a femtosecond x-ray irradiation, the non-adiabatic energy exchange triggers a 
phase transition into low-density liquid phase above the threshold of ~ 0.65 eV per atom in terms of the absorbed dose. 
This semi-metallic state is characterized by a closed band gap, with the local order present in atomic structure. At 
higher doses above ~ 0.9 eV/atom, silicon melts into high-density liquid phase with amorphous atomic arrangement. 
The modeled phase transition occurs within ~ 300 — 500 fs, in a good agreement with the timescales observed in 
experiments. We have also predicted the damage threshold fluence in silicon as a function of the incoming photon 
energy. 

The transition into high-density liquid phase proceeds as a result of the interplay between nonthermal and thermal 
effects. Weakening of interatomic bonds and heating of lattice by excited electronic subsystem triggers ultrafast 
amorphization of silicon. Neglecting electron-phonon coupling results in a significant overestimation of the phase 
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FIG. 10: Damage threshold fluences for silicon corresponding to the low-density liquid and high-density liquid formation as a 
function of photon energy. 



FIG. 11: Nonthermal phase transition: snapshots of atomic positions in silicon irradiated with 10 fs laser pulse of h lo= 1 keV 
photon energy at the absorbed dose of 2.5 eV/atom: (a) t = 0 fs, (b) t = 300 fs, (c) t=0.5 ps, and (d) t = 1 ps. X, Y, and Z 
axes are shown (left-bottom of each panel). 


transition threshold, which then is ~ 2.1 eV/atom. This threshold discrepancy indicates that the non-adiabatic effect 
has to be taken into account in the description of the transitions within x-ray excited silicon. Future experiments with 
FEL should be able to verify these predictions and unveil more details on the radiation-triggered structural transitions 
in silicon. 
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Appendix: Non-adiabatic electron-atom energy exchange 


Boltzmann electron-atom collision integral, I e -at > can be generally written in the following form (similar to (~50|): 


N ~ N poo poo 

E T iJ at = tE / I M e-at(Ei, Ej) I 2 F id (Ei,Ej,E at , E f a ?)dE at dE f a r, 
3=1 i=l J ° 70 


(A6) 


where 


F — 
1 h3 ~ 


{fe{Ei){2 - f e (Ej))f at (E at ) - f e (Ej)(2 - UE^fa^ElrmElr - E at + E, - Ej) ,for * > j 
(fe(Ei)(2 - fe(Ej))f at {E S a ?) - f e {Ej){2 - f e (E z ))f at (E at ))5(E f J n - E at + E j - E t ) ,for i< j, 


is split up into two intervals depending on indices i and j , and summations are running through all the N energy 
states. These state in our case are obtained by diagonalization of the tight binding Hamiltonian. M e _ ot is the 
electron-atom scattering matrix element; f e (E) is the electron distribution function, assumed to be here the Fermi 
function, normalized to 2 accounting for the electron spin; f at is the atom distribution function, taken at the initial, 
E at , and final, states. 

After the integration over the energy conservation condition gives the following: 

poo 

/ F ilJ (E i ,E j ,E at ,E I J n )dE^= (A7) 

Jo 

I f e (Ei)(2 - f e (E 3 ))f at (E at ) - f e (Ej)(2 - f e (Ei))f at (E at - - Ej)), for i > j 

\ f e (Ej)(2 - fe(Ei))f a t(E at - (£, - Ei))0(E at - (Ej - Ei)) - fe(E t )(2 - f e (E 3 ))f at (E at ), for i < j, 


where the 9 -functions are introduced (6(x) = 0 for x < 0, and 6(x) = 1 otherwise). These 6 -functions are allowing 
only for the transitions that conserve energy: in case when an electron jumps up to the final energy level above the 
initial one, j > i, the transitions are only possible as long as the atomic energy, E at , is sufficient to contribute to such 
an event. 

For numerical evaluation of the collision integral over E at fEq. (IA6l) l we assume that the atomic distribution is the 
Maxwellian one with the transient temperature calculated as a kinetic temperature of atoms in a box with periodic 
boundaries: 


/«,(£«,) = 2^. -E. exp (-f£) , 

_ 2 E k in 

at = m at - e ’ 

where E^in is the total kinetic energy of atoms in the super-cell; N at is the number of atoms; T at is the atomic 
temperature in energy units. This assumption is justified as the atomic heating is a slow process compared to 
nonthermal melting, which does not bring the system far out of thermal equilibrium. The initially thermalized atomic 
system only gains kinetic energy, thus increasing the temperature, but it remains in the equilibrium state described 
by the Maxwellian distribution. 

Although one could, in principle, obtain the transient nonequilibrium distribution of atoms by sorting the atomic 
energies into energy intervals, such an approach shows itself as numerically challenging, because an integral and a 
double summation have to be performed in Ea. (IA6l) . In contrast, Maxwellian distribution function for atoms allows 
for exact analytical evaluation of the inner integral. This also significantly speeds up the calculations. 

An integral of the Maxwellian function with the 0-function is as follows: 

gat(E) = fat(Eat)dE at = 2^j■ exp - ( er f _1 ) ’ ( A1 °) 


(A8) 

(A9) 


where erf(x) is the error-function. Note that an integral with the lower limit E = 0 yields 1. 
Combining these equations together, the total collision integral can be written as: 


N 


27r 


N 


E^r^irE We-ati^Ej)? 


3=1 


3=1 


fe(Ei)( 2 - fe(Ej)) - f e (Ej)( 2 - f e (i? l ))Sa t (E i - Ej) , for i > j , 
fe(Ei)( 2 - f e (Ej))g a t(Ej - Ei) - f e (Ej)( 2 - f e (Ei)) ,for i < j , 


(All) 
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FIG. 12: Convergence study of the electron-atom energy exchange rate with respect to the number of atoms in the simulation 
box. Initial conditions applied are: the electron temperature T e = 10000 K, the atomic temperature T a t = 300 K. Constant- 
volume simulations are used here. 


where the integrated Maxwellian function, g at (E), is defined by Ea. (|A10|) . 

The last remaining term to be defined is the electron-atom scattering matrix element M e _ at (E i: Ej). To derive 
it, we use an ab-initio approach proposed by Tully [3(1 37]. Probabilities of non-adiabatic transitions of electrons 
between the energy levels i and j induced by the atomic motion during the current time step, corresponding to the 
atomic displacement R = i?o + SR, can be written in the diabatic representation as follows [59} : 

M e _ at {Ei,E 0 ) = (i(R 0 (t-St))\H e . at (R(t))\j(R 0 (t-6t))), (A12) 

under an assumption of an infinitesimal time-step St, the Hamiltonian can be expanded as H e - at (R(t)) = H e - at (Ro(t— 
St)) + X7H e - at (R 0 (t - St)) • SR: 

M e — a t{Ei, Ej) = {i(R 0 (t-St))\VH e _ at (R 0 (t-8t))\j(R 0 (t-St))) -SR, (A13) 

Utilizing the Hellmann-Feynman theorem, one can rewrite the expression for the potential in terms of the en¬ 
ergy levels (eigenvalues of the Hamiltonian) and the derivatives of the wave-functions with respect to the nuclear 
coordinates: 


(i\VH e - at \j) = (Ej - Ei)(i\V\j) . (A14) 

Thus, for the infinitesimal time-step St, the matrix element can be expressed in terms of the non-adiabatic coupling 
vector (dij = (i|V|j)) |37 : 

<»|V|j>. SR = Sta ■ <?,, = « ((i(t ~ W (A15) 

where R is the atomic velocity; and the wave-functions are defined at the current and the previous step of the 
simulation [35] . The final expression for the matrix element of the non-adiabatic electron coupling to the atoms can 
now be written as: 

M e - a t(Ei, Ej) = ^ ((i(t - 8t)\j(t)) - (j(t - St)\i(t))) (Ej - E t ) (A16) 

where an energy level is taken as a mean value on the current and previous time-steps: Ei = (Ei (t)+Ei (t—St))/2. In our 
case H e _ at (R(t)) is equal to the tight-binding Hamiltonian, thus, the energy levels are the eigenvectors corresponding 
to the eigenfunctions of the tight binding Hamiltonian. These equations (IA16I) and (IA11I) constitute the non-adiabatic 
coupling between the electronic and atomic subsystems which has been used in the current work. 
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The finite-size effects must be investigated prior to any application of the model to a realistic situation. For this 
purpose, we analyzed an electron-atom energy exchange rates in a nonequilibrium model system at the following 
initial conditions: the electron temperature of T e = 10000 K, and the atomic temperature of T at = 300 K. The results 
are shown in Fig. [T2] The figure shows that the number of atoms only slightly affects the electronic temperature, 
and has almost no influence on the atomic temperature. Slight decrease of the electron-atom energy exchange rate 
for low numbers of atoms can be attributed to the contribution of long-wavelength phonons, which appear only for 
sufficiently large simulation boxes. Their contribution is, however, only minor: any differences practically vanish for 
the number of atoms exceeding 216. 
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